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^ ■ ABSTRACT 

\ We present SAUR0N integral field spectrography of the central 1.5 kpc of the nearby 

I/"") , Seyfert 2 galaxy NGC 1068, encompassing the well-known near-infrared inner bar 

observed in the K band. We have successively disentangled the respective contribu- 
tions of the ionized gas and stars, thus deriving their two-dimensional distribution and 
l/"~) ■ kinematics. The [O in] and H/3 emission lines exhibit very different spatial distribution 

and kinematics, the latter following inner spiral arms with clumps associated with star 
^ c"| formation. Strong inwards streaming motions are observed in both the H/3 and [O III] 

fS kinematics. The stellar kinematics also exhibit clear signatures of a non-axisymmetric 

tumbling potential, with a twist in both the velocity and Gauss-Hermite /13 fields. We 
re-examined the long-slit data of Shapiro et al. (2003) using a penalized pixel fitting 
routine: a strong decoupling of the Gauss-Hermite term /13 is revealed, and the central 
decrease of Gauss-Hermite term /14 hinted in the SAUR0N data is confirmed. These data 
also suggest that NGC 1068 is a good candidate for a so-called c-drop. We confirm 
the possible presence of two separate pattern speeds applying the Tremaine- Weinberg 
method to the Fabry-Perot Ha map (Bland-Hawthorn et al. , 1991). We also examine 
jjj ■ the stellar kinematics of bars formed in N-body + SPH simulations built from ax- 

isymmetric initial conditions approximating the luminosity distribution of NGC 1068. 
The resulting velocity, dispersion, and higher order Gauss-Hermite moments success- 
fully reproduce a number of properties observed in the two-dimensional kinematics of 
NGC 1068, and the long-slit data, showing that the kinematic signature of the NIR bar 
is imprinted in the stellar kinematics. The remaining differences between the models 
and the observed properties are likely mostly due to the exclusion of star formation 
and the lack of the primary large-scale oval/bar in the simulations. These models nev- 
ertheless suggest that the inner bar could drive a significant amount of gas down to 
a scale of ~ 300 pc. This would be consistent with the interpretation of the er-drop 
in NGC 1068 being the result of central gas accretion followed by an episode of star 
formation. 

Key words: galaxies: evolution galaxies: individual: NGC1068 galaxies: kinematics 
and dynamics galaxies: Seyfert galaxies: nuclei 
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1 INTRODUCTION 

The fueling of Active Galactic Nuclei (AGN) poses the prob- 
lem of bringing gas in the close neighbourhood of the pu- 
tative central dark mass, a supermassive black hole. Before 
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reaching scales of a few parsecs where turbulent viscosity 
becomes important (Wada & Norman 2002), the angular 
momentum of the gas must decrease by orders of magni- 
tude. QSO are usually associated with a major merging 
event which provides the necessary time varying potential to 
allow this to happen efficiently (Canalizo & Stockton 2001). 
However, the presence of a Seyfert nucleus does not corre- 
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late significantly with the presence of companions or other 
external environmental properties (Maia et al. 2003). In this 
context, departures from axisymmetry in the gravitational 
potential have been advocated to play an important contri- 
bution in the removal of angular momentum of the dissipa- 
tive component (Heller & Shlosman 1994). 

Quadrupole perturbations such as bars are ubiquitous 
in disk galaxies. And indeed bars are good at redistributing 
the gas, and more specifically at concentrating gas within 
their inner regions (Sakamoto et al. 1999). However a cor- 
relation between the presence of a bar and the activity of 
the nuclear region is weak (Laine et al. 2002, Knapen et 
al. 2000, Moles et al. 1995, Mulchaey & Regan 1997, Malkan 
et al. 1998). This is not too surprising since the scales in- 
volved are very different: from the kiloparsec scale bars to 
the presumed central accretion disk. The bar-driven loss of 
angular momentum mainly occurs within the Corotation Ra- 
dius (CR), leading the gas towards the inner resonances, 
e.g. the Inner Lindblad Resonance (ILR) if present, where 
the torques are cancelled out. Gas accumulates at this ra- 
dius, forming an inner ring in which clumps of vigorous 
star formation are often observed (e.g. Schwarz 1981, 1984). 
The next step in moving gas to smaller radii is still under 
great debate. A number of processes, including secondary in- 
ner bars, inner spirals, lopsidedness and minor mergers (see 
Combes 2003 for a review) have been invoked, but none ap- 
pear to provide a necessary and sufficient condition for the 
triggering of nuclear activity. 

Although gas is a very sensitive tracer of the presence 
of a barred potential (e.g. Mundell & Shone 1999), its non- 
linear reponse to even weak non-axisymmetries means that 
it cannot be used directly to derive the gravitational poten- 
tial of the galaxy. In addition, gas flows close to the AGN 
may be dominated by non-gravitational forces due to jets 
and outflow winds (Nelson & Whittle 1996). Therefore, the 
stellar kinematics, although challenging to measure in ac- 
tive galaxies, offer a more direct probe of the underlying 
potential. 

The difficulty in making reliable measurements of the 
stellar kinematics in AGN, particularly in nuclear regions 
complicated by the presence of strong line emission, has re- 
sulted in only a limited number of moderate-resolution, stel- 
lar absorption line studies. Two-dimensional stellar kinemat- 
ics have been published for only a small number of Seyferts 
(e.g. Garo'a-Lorenzo et al. 1997; Arribas et al. 1997; Arribas 
et al. 1999; Ferruit et al. 2004) and in some cases suffer from 
a too restricted field of view for accurate determination of 
the galaxy potential. For most Seyferts, only long slit stud- 
ies (Perez et al. 2000; Emsellem et al. 2001; Filippenko & 
Ho 2003; Marquez et al. 2003) or central velocity dispersion 
measurements (Nelson & Whittle 1995, 1996) are available. 

As one of the closest and most famous Seyfert 2 galax- 
ies, NGC 1068 has been studied at most wavelengths and 
nuclear gas kinematics have been used to constrain its 
overall mass distribution. A good summarizing sketch of 
the observed components can be found in Schinnerer et 
al. (2000, their Fig. 6) where the outer disk and oval, the 
two-arm spiral and the NIR bar are represented. The outer 
oval structure (with a diameter of 90") has been inter- 
preted by Schinnerer et al. (2000) as a primary bar with 
flp ~ 35 km s" 1 kpc" 1 , consistent with the HI ring be- 
ing at its Outer Lindblad Resonance (OLR), and the in- 



ner spiral arms seen in CO corresponding to its Inner Lind- 
blad Resonance (ILR). The NIR bar, which extends up to 
a radius of ~ 16" would then be a secondary bar with 
an estimated pattern speed of Q.^ ~ 140 km s" 1 kpc -1 . 
This is consistent with the recent lower limit estimate of 
135 ± 42 km s" 1 kpc" 1 derived by Rand & Wallin (2004) 
based on an application of the Tremaine- Weinberg (1984) 
method using CO observations of the molecular gas. The 
same authors argued for a different slower pattern speed 
for the CO spiral arms with fip ~ 72 km s" 1 kpc" 1 . How- 
ever, based on the openess of the spiral arms, Yuan & Kuo 
(1998) argued these to be associated with the OLR of the 
inner bar. The only previous two-dimensional stellar kine- 
matics study was achieved by Garcia- Lorenzo et al. (1997) 
who obtained INTEGRAL/WHT spectroscopy in the optical of 
the central 24"x20" of NGC 1068. These authors suggest 
the presence of two puzzling kinematically decoupled stellar 
components in the central 10", offset by about 2'.' 5 from each 
other. Finally, Shapiro et al. (2003) recently constrained the 
velocity dispersion ellipsoid of the disk using high-quality 
long-slit kinematics along the major and minor axes of the 
galaxy. 

In this paper we wish to further constrain the gravita- 
tional potential of NGC 1068 by mapping both its stellar and 
gas kinematics. We have thus used the integral field spectro- 
graph SAUR0N to probe the region of the inner (secondary) 
bar. We report the results obtained from this dataset, anal- 
ysed with the help of N-body + SPH simulations. We first 
describe the SAUR0N observations and associated data reduc- 
tion in Sect. [5] We then briefly present additional data we 
have gathered from different authors (Sect.[3J, and present 
the results in Sect. 2] Numerical simulations for NGC 1068 
are described in Sect. and the resulting kinematics are 
then compared with our data. The results are discussed fur- 
ther in Sect.||j] and summarized in Sect.Q Throughout this 
paper, we will use a distance of 14.4 Mpc for NGC 1068 
(Bland-Hawthorn et al. 1997), leading to a scale of 69.8 pc 
per arcsecond. 



2 SAURON OBSERVATIONS AND DATA 
REDUCTION 

2.1 The SAURON datacubes 

We have observed NGC 1068 with the Low Resolution 
(LR) mode of the integral field spectrograph SAURON during 
a run in January 2002. SAURON delivers about 1500 spec- 
tra simultaneously, homogeneously covering a field of view 
of 41" x 33" with a squared sampling of 0"94 per spa- 
tial element (lens). The spectral domain and resolution are 
[4820 - 5280 A] and 108 km s" 1 (ov) respectively. More de- 
tails on the SAURON spectrograph can be found in Bacon 
et al. (2001), de Zeeuw et al. (2002), and Emsellem et 
al. (2004). The emission lines in NGC 1068 (particularly the 
[O in] lines) are so bright in the central part of NGC 1068 
that the CCD saturates after about 8 min of exposure with 
the SAURON LR mode. We have therefore obtained a 5 min ex- 
posure to probe the central few arcseconds, and three more 
exposures of 30 min to reach sufficient a signal-to-noise ratio 
at the edge of the field of view. 
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Figure 1. Isophotes (with steps of 0.5 magnitude) of the cen- 
tral region as observed by SAURON (stellar continuum, solid lines) 
and with the WFPC2 camera aboard HST (F547M filter, dashed 
lines), after (seeing-)convolution by a gaussian of 2". The agree- 
ment is excellent except at the border of the HST/PC field where 
edge effects due to the convolution are seen. 



2.2 Data reduction with XSAURON 

The data reduction of all four exposures was achieved us- 
ing the dedicated XSAURON software, and an automated 
pipeline available within the SAURON consortium (see Ba- 
con et al. 2001, and de Zeeuw et al. 2002 for details). The 
main steps include: removal of the CCD signature, extrac- 
tion of the spectra using a mask built from an optical model 
of the telescope+spectrograph, wavelength calibration, low 
frequency spectral fielding, cosmic rays removal, homoge- 
nization of the spectral resolution in the field, subtraction 
of the sky contribution using 146 dedicated lenses (1.9' away 
from the central field), and flux calibration. The sky spec- 
tra were carefully checked for any significant contribution 
(emission/absorption) from NGC 1068 itself. 

The flux-calibrated individual exposures are then accu- 
rately centred relatively to each other using the associated 
reconstructed images, and merged. Since the exposures are 
slightly offset from each other, we use the dithered datacubes 
to spatially resample the datacubes to 0'.'8 per spaxel (spa- 
tial pixel). Before merging, the spectra which exhibited sat- 
urated pixels were removed from the datacubes. This prac- 
tically means that the spectra in the central 4" of the final 
merged datacube are only coming from the 5 minutes un- 
saturated exposure. At the edge of the SAURON field of view, 
the signal-to-noise ratio gets down to about 30 per pixel, 
sufficient for our goal of deriving the stellar kinematics. 



2.3 Coordinate centering 

We used the F547M WPFC2/HST (extracted from the 
ESO/ST-ECF archive) data to determine the relative spa- 
tial position of our SAURON datacube (see Fig. . Capetti 



et al. (1997) provided the absolute coordinate of the cen- 
tral peak of this optical image (a — 02h42m40.s711, S — 
-00deg00'47."81 - J2000, FK5, 80 mas accuracy), to be com- 
pared with the putative position of the central engine (CE) 
identified as the SI radio source by Muxlow et al. (1997) 
at a = 02h42m40.s7098, S = -00deg00'47."938. Taking SI 
as our reference (0,0), we thus aligned the SAURON recon- 
structed image with the F547M WFPC2 exposure where 
the peak was assumed to be 0'.'13 North and 0'.'02 East of 
the CE (including a potential rotation of the SAURON field of 
view). This ensures an absolute positioning of our SAURON 
datacube on SI with an error better than O'.'l in translation 
and than 0.5° in rotation (the latter being dominated by the 
uncertainty on the relative angle between the WFPC2 and 
SAURON data). All maps presented in this paper are oriented 
in the classical way, with North up, and East left. 

2.4 Stellar kinematics 

The stellar contribution to the SAURON spectra is highly 
contaminated by strong emission lines: H/3A4861, the 
[Om]A4959,A5007 and [N i]A5198,A5200 doublets. We first 
identified the spectral regions which are significantly con- 
taminated by emission. This required a detailed examination 
of the spectra in the merged datacube, particularly in the 
central 5" where the emission lines are very strong and wide. 
A first estimate of the stellar kinematics is then derived ex- 
cluding the contaminated pixels. We achieve this by using 
a direct penalized pixel fitting routine (pPXF, Cappellari 
& Emsellem 2004): the algorithm finds the mean velocity 
V and velocity dispersion a which minimizes the difference 
between the observed galaxy spectrum and the spectrum of 
a stellar template convolved by the corresponding gaussian 
(of mean V and root mean square a). This is performed 
with the galaxy and template star spectra rebinned in In A. 
We use a single template star for this first estimate, namely 
the K2 giant HD 26162, also observed with SAURON. Best-fit 
values for V and a are thus obtained at each position (for 
each lens/spectrum) independently. 

We then use this initial estimate of the stellar kine- 
matics to derive an optimal stellar template for each in- 
dividual spectrum. This is usually done using a large li- 
brary of stellar templates. In the case of NGC 1068, and be- 
cause of the significant number of pixels we had to exclude 
from the fit, this would make the fitting process strongly 
degenerate. We therefore decided to restrict our stellar li- 
brary to include three different stellar templates from the 
single-age single-metallicity stellar population (SSP) models 
of Vazdekis (1999); this proved sufficient to properly fit the 
underlying stellar contribution (see Fig.^. A multiplicative 
polynomial with a maximum degree of 6 was included in the 
fit in order to account for small residual differences between 
the stellar libraries and the SAURON spectra. 

We finally iterate by measuring the velocity V and dis- 
persion a, as well as the third and fourth Gauss-Hermite 
moments /13 and hi with the pPXF routine, this time us- 
ing the optimal templates obtained from the previous step. 
/13 and hi are indicators of the skewness and peakiness of 
the line-of-sight velocity distribution (relatively to a Gaus- 
sian). Although the signal-to-noise ratio per pixel (> 30) is 
sufficient everywhere in the SAURON field of view to derive re- 
liable /13 and hi values, these should be taken with caution 
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Figure 2. Fits of the SAURON spectra of NGC 1068 using the Vazdekis stellar library. Six spectra (black lines) and their corresponding 
fits (red lines) at different locations in the SAURON field are shown, the individual spatial locations being indicated by arrows on the 
reconstructed intensity SAURON map. The vertical position of the best-fit spectra are arbitrary shifted for legibility, the residual spectra 
being presented in the same panels (blue lines). The SAURON map has North up, and East left. 



considering the limited spectral domain of the SAURON dat- 
acube, the medium spectral resolution of the SAURON spectra 
(a = 108 km s" 1 ) and more importantly the contamination 
by bright emission lines. A comparison with long-slit data is 
helpful in this context (see Sect. 14.21 . 



2.5 Gas distribution and kinematics 

The spectra resulting from the fitting procedure described in 
the previous Section were subtracted from the original data 
to provide pure emission lines spectra. The wavelength range 
of our observations includes the H/3, [O iii]AA4959,5007 and 
[N i] AA5198,5200 emission lines. These five lines are detected 
over our complete field of view, although the lines of the [N i] 
doublet were significantly weaker than the H/3 and [O ill] 
lines and were barely detected in off-nuclear regions. The 
basic parameters of these emission lines (intensity, centroid 
velocity and velocity dispersion) were derived from Gaussian 
fitting of their profile using the fit/spec software (Rousset 



1992). When relevant, the lines were constrained to have 
fixed or bounded ratios ([Om]A5007 / [Oni]A4959 = 2.88; 
0.7 < [Nl]A5198 / [Nl]A5200 < 2.0) 1 . To stabilise the fit 
of the weak [N i] lines, we forced them to share the same 
velocity and width as the H/3 line. 

Careful examination of the spectra showed the presence 
of several kinematically distinct components in the profiles 
of the H/3 and [O ill] emission lines, in agreement with the 
results of Arribas et al. (1996) and Garcia- Lorenzo et al. 
(1999). In our data, we have identified 3 different systems: 

• a first component, which is observed everywhere in our 
field of view. It is relatively narrow (typical dispersion of 
100 km s~ x or lower beyond 5" from the nucleus) everywhere 
except in the nuclear regions, and it is therefore termed "nar- 
row" hereafter. Despite its broadening in the nuclear region, 

1 Range estimated using the Mappings Ic software with an 
electronic temperature of 10 4 K and electronic densities from 
0.1 cm -3 to 1000 cm -3 (Ferruit et al. 1997). 
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Figure 3. Examples of spectra of NGC 1068 obtained with SAURQN. The spectra, normalized to a peak intensity of unity, have been 
truncated to show only the wavelength range including the H/3 and [O III] AA4959, 5007 lines. The dotted vertical lines indicate the expected 
position of a line at rest in the galaxy frame (systemic velocity of 1144 km s —1 ). (1) Nuclear spectrum. (2) and (3) Spectra located at 
positions (— l"3,+l'/2) and (+0'.'5,— 2'.'8) where both the narrow and additional components are present. In panel 2 the narrow component 
is redshifted with respect to the systemic velocity (the additional component being blueshifted), whereas in panel 3 it is the opposite. (4) 
Spectrum located at (+5'.'3,+8'.'4) in which both the narrow and broad components are present in the [Om] lines. (5) and (6) Spectra 
located (+13'.'3,— 6'.'8) and (— 16'.'3,+14'.'8), respectively and in which only the narrow component is observed (but with very different 
[Om]A4959 / H/3 ratios). 



and the presence of the other systems, it is possible to follow 
this "narrow" component down to about 2" from the centre. 

• a broad (dispersion higher than 600 km s" 1 ), 
blueshifted component (hereafter termed "broad"), which 
is observed up to 8"-10" from the nucleus and often appears 
as a blue wing in the [O ill] line profile (see e.g., panel 4 of 

Fig. B) . 

• a very spatially localized, intermediate-dispersion com- 
ponent (hereafter termed "additional") observed in the 
vicinity of the nucleus and either strongly blueshifted 
(North-East of the nucleus) or redshifted (South- West of 
the nucleus, see e.g., panels 2 and 3 of Fig. £0 with respect 
to the systemic velocity of the galaxy. 

Examples of spectra displaying these three kinematic 
components are shown in Fig. E3 which can be compared 
with Fig. 7 in Garcfa-Lorenzo et al. (1999). Our narrow, 
broad and additional components correspond to components 
(1), (2) + (3) and (4a) + (4b), respectively, in Garcfa-Lorenzo 
et al. (1999). In our data, it was not possible to disentangle 
their components 2 and 3. This is probably due to our spatial 
sampling (0'.'8 per lens), which makes it difficult to study the 
(complex) central few arcseconds where component 3 is iden- 
tified by these authors. For the same reasons, conducting a 
similar comparison with the decomposition used by Arribas 
et al. (1996) proved very difficult because their observations 



cover only the nuclear regions (see their Sect. 3.2.1). It must 
also be emphasized that each system identified in the SAUR0N 
spectra may itself result from a blend of separate (and unre- 
solved) velocity systems: this is for instance clearly the case 
in the central few arcseconds (Cecil et al. 2002). 

At distances > 5" from the nucleus (corresponding to 
1822 spectra out of a total of 1943), the line profiles were 
simple enough for an automated fit to be conducted. The 
various maps inferred from the results of this automated 
fit were carefully checked and only a small number of spec- 
tra had to be fitted a second time individually. In contrast, 
the automated fit to spectra in the inner regions (radii less 
than 5") was unstable, as expected. All spectra had to be 
examined and fitted individually and quite often additional 
constraints (especially on the width of the broad and addi- 
tional components) were imposed. Our results in these inner 
regions (which, however, represent only a small fraction of 
our field of view) are therefore less reliable than those for 
the outer regions. 

Given the known kinematic complexity of the nuclear 
regions (Cecil, Bland & Tully 1990) and the limitations of 
our dataset in this region, we do not discuss the properties 
of the additional and broad components further. In the fol- 
lowing, we focus on the properties of the narrow component, 
which is more representative of the underlying galaxy. 
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Figure 4. Schematic of the structures observed in NGC 1068. Left panel: 10'xlO' DSS I band image. The central 40" are saturated. 
Isophotes are shown with a step of 0.5 mag.arcsec 2 . Middle panel: sketch at the same scale as the DSS image. The major-axes of the outer 
disk, outer oval and inner near- infrared bar are indicated (dashed lines). The HI kinematic axis is shown as a dotted line for comparison. 
Right panel: zoomed sketch (xlO) showing the location of the near-infrared bar and the CO arms (see Schinnerer et al. 2000, Fig. 1 and 
6). The extent of the SAUR0N field of view is indicated by a dashed polygon. 



Errors on the emission line parameters were derived by 
use of Monte Carlo realizations repeating the fitting proce- 
dure 500 times using simulated emission line spectra. These 
were built as the sum of a synthetic noise-free spectrum 
(with spectral characteristics typical of the observed spec- 
tra) and noise. Additional errors on the centroid velocities 
(18 km s" 1 for 3a calibration errors) and the Full Width at 
Half Maximum (FWHM, 0.1 A peak to valley for variations 
of the "instrumental" FWHM over the field of view) were 
included in the final error budget. In the case of the centroid 
velocities of the narrow component, the contribution of the 
fit to the total error was negligible for spectra dominated 
by this component (i.e. typically at radii > 5") and with 
a peak signal to noise ratio larger than 10 (i.e. over most 
of our field of view): the overall accuracy on the measured 
centroid velocity is then better than 20 km s _1 . 



3 ANCILLARY DATASETS 

In this Section, we briefly describe ancillary datasets we use 
in this work, namely some optical and near-infrared images 
as well as the HI velocity rotation curve for the dynamical 
modelling (Sect.[SJ , Ha, and CO distribution and kinematics 
for comparison with our emission line SAUR0N maps. 



3.1 Ground-based photometry 

A deep B band image was used to probe the outer disk of 
NGC 1068 up to a radius of 200". This image was obtained 
with the lm telescope on Mount Laguna (Cheng et al. 1997) 
and is a combination of 3 exposures of 300s each. It has a 
field of view of about 800" sampled at tf.'A per pixel. We 
also made use of the Digitized Sky Survey I band (available 
via the ESO/ST-ECF archive: http://archive.eso.org/) and 
2MASS K band images, both having a scale of about 1" per 
pixel. Finally, a high resolution K band image was obtained 
by Peletier et al. (1999): this has a pixel size of 0'.'248, a field 
of view of about l' 2 , and a seeing of 0'.'5 (FWHM). 



3.2 HI velocity curve 

The HI velocity curve we used has been published in 
the Ringberg Standards (RS hereafter, Bland-Hawthorn et 
al. 1997) from the work of Brinks et al. (1997). It comprises 
measurements up from the centre to about 200" , the last ra- 
dius at which Brinks et al. (1997) detected the low surface 
brightness HI disk. The HI rotation curve decreases from 
about 130 km s" 1 at 30" to 95 km s^ 1 at 180" from the 
centre. 



3.3 The Ha distribution and kinematics 

Dr. J. Bland-Hawthorn provided us with the Ha Fabry Perot 
data as published in Bland-Hawthorn et al. (1991), including 
the luminosity and mean velocity maps. The datacube was 
obtained with the Hawaii Fabry-Perot Interferometer (HIFI) 
at the CFHT (Mauna Kea). The velocity resolution was 
65 km s , the spatial sampling and resolution (FWHM) 
being 0'.'43 and 0'.'8 respectively. The Ha velocity field was 
published in Dehnen et al. (1997). 

3.4 The CO distribution and kinematics 

We also make use of the high resolution CO interferometric 
data published in Schinnerer et al. (2000, S+00 hereafter). 
This 12 CO (1-0) dataset has been obtained with the IRAM 
millimeter interferometer on the Plateau de Bure (France), 
providing a circular beam of l'.'4 (sampling of 0"4), with a 
channel width of 10 km s _1 . Apart from a central ring-like 
structure, the molecular gas is mainly distributed along a 
two-arm spiral just outside the inner near-infrared bar, and 
which is associated with the ILR of the outer oval (S+00). 

3.5 Long-slit stellar absorption kinematics 

Finally, we include results from long-slit spectroscopy con- 
ducted by Shapiro et al. (2003) who kindly made the fully 
reduced datasets available to us: stellar kinematic profiles 
along the major and minor-axes of the galaxy (PAs of 80 
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Figure 5. Reconstructed maps of the stellar kinematics of NGC 1068 obtained from the SAURQN data cube. From left to right: Intensity, 
velocity V, velocity dispersion a (top panels), and higher order moments hz, and /14 (bottom panels). The thick contour approximately 
marks where the kinematic measurements are very significantly perturbed by the presence of very strong and wide emission lines. Boxes 
in the top right corner of the panels show the ranges covered by the color bars. 



and 170 degrees) were obtained with a slit width of 3" cen- 
tred on the Mgb triplet around 5175A. The stellar kine- 
matics published in Shapiro et al. (2003, hereafter Sh+03) 
were obtained using a Fourier fitting algorithm. In this pa- 
per, we reanalysed these data favouring the pPXF routine as 
it allows an optimal selection of regions uncontaminated by 
emission lines, and to derive robust estimates of moments up 
to hi and /14. We have only included wavelengths between 
5237 and 5549A to avoid contamination from the very bright 
emission lines present in these spectra of NGC 1068. 



4 DATA ANALYSIS 

4.1 Morphology and position angles 

In this Section, we derive new values for the main morpho- 
logical parameters of the different components of NGC 1068 
(bar, oval, outer disk), since these are important parameters 
when comparing different datasets and models. A schematic 
summarizing the structures observed in NGC 1068 is pro- 
vided in Fig. 0] 

We first wish to reexamine the position angle (PA) and 



ellipticity of the outer disk. The RS quote a value of 106° ±5, 
which in fact corresponds to the average position angle of the 
kinematic axis as fitted on the large scale HI data (Brinks 
et al. 1997). The total HI surface brightness in the outer 
disk (100-200") is rather low and exhibits a North/South 
asymmetry at its outer edge (lower surface brightness in 
the South). At this radius the stellar component is easily 
observed using the DSS I band image: outside 190", there 
is a clear mildly flattened component with a PA between 
74.5 and 84°, therefore at least 20° from the kinematic PA 
measured in HI. A value of 80° ± 5 is consistent with the op- 
tical and smoothed HI surface brightness shown in Dehnen 
et al. (1997), and we will adopt this value as the apparent 
photometric PA of the outer disk, thus different from the 
HI kinematic PA. The axis ratio of the best fit ellipse of the 
optical outer disk is between 0.8 and 0.85. 

We then remeasured the characteristics of the outer oval 
and near-infrared bar using ellipse-fitting. Both the DSS I 
band and the B deep images lead to a radius of 90", an axis 
ratio of 0.8 and a PA of 5° for the outer oval, perfectly con- 
sistent with the RS values. Using our high resolution K band 
image, we find an average position angle of 44.5° ± 0.5° for 
the near- infrared bar (between 10 and 16" radius), as com- 
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pared to the RS value of 48° ± 2° (from Scoville et al. 1988). 
Our value is however consistent with the one provided by 
Thronson et al. (1989) of 45.0° ± 0.5°. In the following, we 
will therefore use an average value of 44.5°. The minimum 
axis ratio is 0.45 at a radius of 15'.'5. 



4.2 Stellar kinematics 

The SAURON stellar kinematics are presented in Fig. [S] Inside 
the central 4", the measurements are significantly perturbed 
due to the emission line contamination and should be taken 
only as indicative. The stellar velocity field clearly exhibits 
strong departures from axisymmetry, with an "S" shaped 
zero velocity curve, and the line of maximum velocity hav- 
ing a changing position angle. The amplitude of the velocity 
field reaches ~ 115 km s _1 at a radius of about 10". We 
wish to draw attention to a small perturbation of the order 
of 20km s _1 at the South-East edge of the field (absolute 
values being smaller, around a band going from 12" East, 
6" South to -18" East from the nucleus; see Sect. 14.31 and 
a similar trend is mirrored on the opposite side. The ve- 
locity dispersion map shows a rise towards the center with 
a ~ 60 km s _1 at 20" from the centre along the major-axis 
and between 100 and 200 km s _1 in the central 10". There 
is a slight asymmetry in the dispersion map (dispersions be- 
ing higher by ~ 20 km s _1 on the western side of the field), 
probably the result of a residual gradient of the spectral res- 
olution over the field of view: however, this does not affect 
our conclusions. The ft 3 map displays a significant change 
of sign from the SE quadrant to the NW quadrant. The 
structure of positive and negative /13 is roughly elongated 
along the position angle of the NIR bar. There is an elon- 
gated ring-like region of positive values in the hi map (with 
maxima around 0.06-0.1 between radii of 8 and 12") and 
a depression inside a radius R of 8" with hi going slightly 
negative. We cannot however follow this inward decrease of 
hi for R < 4" because of the emission-line contamination. 

We now compare the kinematics published by Sh+03 
and obtained via a Fourier fitting technique, with our reanal- 
ysis of the same dataset using pPXF (Fig. . The central 
velocity dispersion does not peak so clearly in the profiles 
reextracted from the Sh+03 data: we were unable to repro- 
duce the high central a value even by changing our stellar 
template or spectral domain. There is also a slight flattening 
of the major-axis velocity gradient in the central 2" in the 
Sh+03 profiles, which we do not see in our measurement. 
The hz profiles are anticorrelated with V with peak values 
of ±0.15 around 10" along the major-axis. The hi profiles 
exhibit a significant depression in the central 5" with a nega- 
tive minimum of ~ —0.04 at the centre. The maximum value 
of hi is reached at a radius of about 14" along the major-axis 
and 10" along the minor-axis confirming the elongation of 
this ring-like structure. The presence of a drop in the central 
dispersion profile is confirmed by Shapiro & Gerssen (private 
communication) who reexamined this dataset and applied 
their own pixel fitting routine. We presume that the central 
peak in the stellar velocity dispersion observed in Sh+03 is 
due to the influence of the strong [NI] doublet when using a 
Fourier fitting program to extract the kinematics. Consider- 
ing this reanalysis, NGC 1068 is therefore a new candidate 
for the presence of a so-called a-drop (Emsellem et al. 2001; 
Wozniak et al. 2003; Marquez et al. 2003). 
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Figure 6. Comparison between stellar kinematic profiles ex- 
tracted from the SAURON (solid lines) and Sh+03 (crosses). From 
top to bottom: mean velocity, velocity dispersion, /13 and /14. The 
SAURON data have been averaged over a slit width of 3" 2 (the cen- 
tral 3" are unreliable due to the presence of dominating emission 
lines). The mean velocity and dispersion profiles originally pub- 
lished by Sh+03 are provided as circles (with error bars). The 
"major-axis" here corresponds to a PA= 80°, and the "minor- 
axis" to PA= 170°. 



Overall the SAURON stellar kinematics averaged over a 
reconstructed slit of 3" 2 width compare reasonably well with 
the published long-slit data of Sh+03, although there are 
significant discrepancies worth mentioning: the SAURON dis- 
persion values are too high on the western side of the major- 
axis, confirming the fact that the SAURON data have a resid- 
ual spatial variation of the spectral resolution. In contrast 
to the profiles measured from the Sh+03 data, our /13 data 
do not show a significant gradient. This is primarily due to 
our lower spectral resolution (a = 108 km s - , as compared 
to about 30 km s" 1 for Sh+03), which obviously also affects 
our hi measurements. 

We do not see any hint for the decoupled kinematic 
structure claimed by Garcia-Lorenzo et al. (1999, 1997). 
Although the derived stellar kinematics within the very 
central 4" are unreliable, we should be able to detect the 
presence of a large East- West velocity gradient as implied 
by the pinched isovelocities in the maps of Garcia-Lorenzo 
et al. (1997). There is also no hint of an abrupt veloc- 
ity change in the minor-axis kinematics of Sh+03. This 
strongly suggests that the kinematic structure observed by 
Garcia-Lorenzo et al. (1999, 1997) is an artefact in the 
INTEGRAL / WHT data. 
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Figure 7. SAURQN maps of the gas distribution and kinematics for NGC 1068. Top row (from left to right): distribution of the total, and 
narrow components, and velocity and velocity dispersion for the H/3 narrow component. Bottom row: same as top row, but for the [Om] 
component. Boxes in the top right corner of each panel show the ranges covered by the color bars. 



4.3 Gas distribution and kinematics 

The maps of the distribution and kinematics of the H/3 and 
[O in] components, as reconstructed from the results of mul- 
tiple Gaussian fits to the data, are displayed in Fig. [7| 

H/3 emission from the narrow component is ubiquitous 
in our field of view but the brightest region corresponds to 
the inner 3" of the galaxy. Away from these nuclear regions, 
the narrow-component H/3 emission is dominated by the con- 
tribution from the spiral arms. The distribution of H/3 shows 
the known northern and southern spiral arms although they 
seem almost connected at this resolution. Comparison be- 
tween the total (all kinematic components included) and 
narrow-component-only H/3 maps (top left panels of Fig. [7J 
show little difference between the two maps outside from the 
nuclear regions, outlining the fact that the narrow compo- 
nent dominates the H/3 emission away from the nucleus. 

As for H/3, emission from the inner 3" of the galaxy 
also dominates the narrow-component [O in] map. However, 
away from the nucleus the distribution of the [O ill] emission 
differs significantly from the distribution of the H/3 one. The 
distribution of the [O ill] emission is very asymmetric, thus 
found predominantly North-East of the nucleus, and does 
not trace the spiral arms. Instead, it traces the northern 
ionization cone (see e.g. Pogge 1988). 

Despite the differences between the distribution of H/3 
and [Om], the overall morphologies of the H/3 and [O ill] 
velocity fields (see Fig. [7J are very similar. They both dis- 
play the prominent "S" shaped structure, which is already 
present in the data of Cecil et al. (1990; see also Sect. I4.4jl 
and of Garci'a-Lorenzo et al. 1999 (their Fig. 16). However, 
if we take a closer look at these velocity maps, we see dif- 



ferences between the observed [O ill] and H/3 velocities in 
some regions. We have therefore built a map of v([Om]) 

— v(H/3), which is displayed in Fig. |H1 Differences up to 
±150 km s~ are measured, with uncertainties of typically 
40 km s _1 (3<t). Strong positive values of observed v([0 ill]) 

— v(H/3) are located in the South-East half of our field of 
view, while strong negative values are found on the other 
half. These qualitatively follow the perturbations observed 
in the stellar kinematics (see Sect. I4.2ll . 

We also show the velocity dispersion maps of H/3 and 
[O ill] (right panels of Fig. |7J . In both maps and if we ex- 
clude the nuclear regions, the regions of bright emission dis- 
play a dispersion lower than their lower-surface-brightness 
surroundings. We checked that this was not a bias intro- 
duced by differences in signal-to-noise ratio between bright 
and weak regions 2 . The signal to noise ratio in our H/3 and 
[O in] observations is > 20 over most of our field of view and 
our Monte-Carlo simulations (see Sect. 1231k show that, in 
this case, our relative uncertainties on the measured disper- 
sion values are < 10 % (3a). We observe differences of the 
order of 0.5 A for dispersion values of typically 2.5 A and 
this trend is therefore real. 

Last, we have built the map of the [Oiii]A5007 / H/3 
ratio (see Fig.|5Jl. The range of variation of this ratio is ex- 
tremely large with H/3 dominating in the arms (ratio below 
unity) and [O ill] dominating in the ionization cone (ratio 
reaching values up to 13-14). The peak of high [Om] / H/3 



2 There is a systematic trend to obtain larger dispersion values as 
a result from the fit when the signal to noise ratio in a spectrum 
decreases. 
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Figure 8. Reconstructed map of the difference between the cen- 
troid velocities of the [O III] and H/3 lines of the narrow com- 
ponent, with contours of the H/3 narrow-component map overi- 
mosed. Contours: 10, 20 and 40 erg s — 1 cm -2 . 
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Figure 9. Reconstructed map of the [O in]A5007 / H/3 ratio, with 
contours of the H/3 narrow-component map overimosed. Contours: 
10, 20 and 40 erg s _1 cm -2 . A clipping has been applied to the 
map, putting to pixels where the relative accuracy of the ratio 
was worse than 20% (only a few pixels south-west of the nucleus). 



ratio around a PA of 30° corresponds to the northern ion- 
ization cone. Note that the presence of the well-defined peak 
suggests that, within the cone, the excitation of the gas de- 
creases from its center to its edges. This result must however 
be taken with caution as the regions responsible for this peak 
are located at the very edge of our field of view. 

There does not seem to be a systematic association be- 
tween regions which exhibit large velocity differences be- 
tween [O in] and H/3 (Fig. |HJ , and those with high [O in] / 
H/3 ratio (Fig.[5J. Outside a radius of 5", we observe large ve- 
locity differences up to 100 km s _1 and rather normal [O ill] 
over H/3 line ratios. This suggests that the observed kine- 
matic disturbance is not directly linked with the AGN ac- 
tivity. 

4.4 Comparison with the Ha and CO maps 

The agreement between the SAUR0N H/3 flux and velocity 
maps with the corresponding Ha maps, as shown in Fig. HOI 
is remarkable considering the very different instrumental 
setup, and the fact that dust extinction is present (Bruh- 
weiler et al. 2001). Apart from the central few arcseconds, 
where scattering effects are important and AGN related 
emission is present, extinction significantly affects H/3 (as 
probed by a ratio of Ha over H/3 well over 3; Bruhweiler 
et al. 2001) in the South- West part of the arm (SW clump 
at a radius of about 14"; Bruhweiler et al. 2001). Both the 
Ha and H/3 velocity field exhibit a very strong spiral-like 
perturbation with the zero velocity curve displaying a very 
wavy shape. The Ha map is analysed further in Sect. 14. fl 

A comparison between the maps of the 12 CO (1-0) line 
flux and the corresponding Ha emission is shown in Fig. 1101 
They both roughly follow the two-arm spiral structure with 
a diameter of about 40" (often mentioned as the circumnu- 
clear ring). There are however some very significant differ- 
ences which we emphasize here. Firstly, the Ha emission is 
much more asymmetric with respect to the centre with a 
brighter northern arm. Secondly the CO spiral is clearly off- 
set from the Ha arms: it is on the inner side of the SW Ha 
clump, but seems to lie outside the Ha arm in the North- 
West. A similar comparison at higher spatial resolution us- 
ing the HST/WFPC2 images (Bruhweiler et al. 2001) shows 
that the spiral arms seen in CO and Ha are in fact offset 
from each others (the northern and southern CO arms being 
outside and inside the corresponding Ha arms, respectively). 
This is commonly observed in spiral barred galaxies (Sheth 
et al. 2002). In the case of NGC 1068, the offset could be 
partly explained by the fact that most of the Ha emitting 
regions are associated with starbursting HII regions. Davies 
et al. (1998) thus showed that most of the young stars in 
the ring formed in compact clusters in a relatively recent 
short burst within the last 30 Myr, about ten dynamical 
timescales at the radius of the ring. But extinction also plays 
a role here since Ay of ~ 1 — 3 mag have been measured 
for these clusters (Davies, Sugai & Ward 1998). We there- 
fore do not expect a perfect coincidence between the CO 
and Ha emission line gas (according to Davies et al. 1998, 
the ionization cone cannot have a significant effect on the 
star formation). Since the near side of the disc is South of 
the nucleus, we cannot discard the possibility that we only 
clearly detect HII regions closer to the edge and in front of 
the molecular arms, which would explain the relative loca- 
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Figure 10. Top panels: comparison between the SAURON H/3 flux (contours, left panel) and velocity (contours, right panel) and the 
corresponding Ha maps obtained by Bland-Hawthorn et al. (1991) with a HIFI Fabry-Perot. The velocity step used for the SAURQN 
isovelocities is 25 km s _1 . Bottom panels: comparison of the Ha maps with the CO distribution (contours: Schinnerer et al. 2000). Boxes 
in the top right corner of each panel show the ranges covered by the color bars. 



tions of the ionized and molecular arms. We however do not 
have a good tracer of the currently ongoing star formation 
compared to the slightly older HII regions, and it is therefore 
not possible to distinguish between the scenarios of dissoci- 
ation of molecular gas due to hot stars, or segregation of CO 
and HII regions (as seen in the spiral arms of M83, M100, 
or M51, see e.g. Rand et al. 1999) due to a spiral density 



wave, or even some other unknown effect: observations of 
mid-infrared emission lines are required to test this. 



4.5 Harmonic Analysis of the Ha Velocity Field 

Our SAURON H/3 and [O ill] velocity fields, as well as the Ha 
one, exhibit complex kinematic features such as "S" -shape 
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Figure 11. Left panel: observed Ho velocity field. Top row, from left to right: best fit circular velocity component alone, corresponding 
residual map (Data - V ro t model) and the derived rotation curve assuming an inclination of 40°, PA= 87°, and V S ys= 1140 km s _1 . 
Bottom row, from left to right: harmonic reconstruction including the rotational velocity component, corresponding residual map (Data - 
Harmonic Decomposition Fit) and the derived radial velocity component normalized with v* = V ro t sin(i) as a function of galactocentric 
radius. The solid curves in the bottom right panel show the radial velocity component according to the analytic spiral model as described 
in Sect. l4*5l The best matching spiral model with CR at 31" corresponding to the spiral pattern Q p = 109 ± 5 km s _1 kpc -1 is shown 
in red. Radial velocities for spirals with pattern speeds of 99 and 119 km s~ 1 kpc -1 are shown in green and blue colours respectively. 



or wobbles in the zero velocity curve. We aim to extract 
the effect of the prominent NIR bar and that of the spiral 
arms from the velocity field of NGC 1068. Given that the 
Ha field covers a significantly larger area than the SAURON 
field (at the expense of a shorter spectral coverage), we will 
therefore apply our analysis method on the former. Assum- 
ing that circular motion is the dominant feature, and that 
there is no strong warp in our field, we carry out an analysis 
based on the harmonic decomposition of the line-of-sight ve- 
locity Vios (Schoenmakers, Franx & de Zeeuw 1997; see also 
Franx, van Gorkom & de Zeeuw 1994). This formalism im- 
plies expanding the Vi os field into harmonic series, where the 
first terms (cos 9 and sin 9) are the rotational and radial ve- 
locity components, and the higher harmonic terms (cosmf? 
and sin m9) provide information about perturbations on the 
gravitational potential (Wong, Blitz & Bosma 2004; Fathi 
2004). The first, second, and third harmonic components, 
are sufficient for studying specific elements of the pertur- 
bations on the underlying potential. A perturbation of or- 
der m creates m — 1 and m + 1 line-of-sight velocity terms 
(eg. Canzian 1993 and Schoenmakers, Franx & de Zeeuw 
1997), i.e. third harmonic terms contain information about 
an m = 2 bar or a two-arm spiral perturbation. 

To obtain the kinematic PA, the systemic velocity, and 
circular velocity contribution to the Ha Vios as a function of 
galactocentric radius, we apply a tilted-ring method similar 
to that by Begeman (1987). We first assume an inclination of 
40° , then for each galactocentric radius, we find the best fit 
PA and systemic velocity. These parameters did not show 
a significant variation throughout the field: the PA varia- 
tion is found to be less than 10° and the Kys varies less 



than 20km s 1 . Following the standard procedure, we then 
fix the PA and V sys to the average values. This step yields 
PA= 87° and V sys = 1140 km s _1 , which are consistent with 
the values mentioned in the rest of the present paper. Keep- 
ing these parameters fixed, we now fit the rotational com- 
ponent, and derive the rotation curve presented in Fig. 1111 
Reconstruction of the circular velocity component and sub- 
traction from the observed velocity field yields the residual 
non-circular velocities as presented in Fig. ^3 

We explore the non-circular velocity components by fit- 
ting the residual map with the higher order harmonic terms 
up to and including order 3. The second harmonic terms are 
mainly consistent with zero. The first and third harmonic 
terms show a behaviour similar to that of an analytically 
derived logarithmic two-arm spiral perturbation (Wong et 
al. 2004). We therefore construct a library of velocity fields 
perturbed by two armed spirals with different spiral struc- 
ture characteristics. The library of models include a range of 
pitch angles, perturbation amplitudes, and spiral arm sizes. 
We compare the harmonic terms and non-circular velocity 
features with those of the models to find a model which 
best resembles the observed features. We find that the spi- 
ral model with a pattern speed Q ~ 109 ± 5 km s _1 kpc -1 , 
a pitch angle of 15°, a spiral amplitude of 0.17, and with its 
corotation radius at 31" provides the best fit to the observed 
Ha velocity field. In the region of the inner near-infrared 
bar, there is however significant disagreement between this 
model and the observed maps, the spiral model alone being 
obviously inadequate for describing the observed features. 



4.6 Application of the Tremaine-Weinberg 
method 

Rand & Wallin (2004; RW04 hereafter) recently argued for 
the presence of two different pattern speeds for the NIR in- 
ner bar and the spiral arms which surrounds it in NGC 1068. 
They quote a value of Q p ~ 135 ± 42 km s -1 kpc -1 for the 
inner bar and a lower value of SI® ~ 72 ± 4 km s _1 kpc -1 
for the outer spiral arms via the Tremaine-Weinberg method 
(Tremaine & Weinberg 1984) applied to the CO molecular 
line data. This estimate for the pattern speed of the inner 
bar should be taken as a lower limit, considering the pres- 
ence of a second outer tumbling component. The validity of 
using this technique on a gaseous component was tested by 
RW04 with Nbody + SPH simulations. A similar test was 
also recently performed by Hernandez, Wozniak, Carignan 
et al. (2005) who emphasized the need to avoid regions of 
shocks. We attempted to estimate the pattern speed of the 
spiral and bar structures in NGC 1068 via the same tech- 
nique but on the large-scale Ha maps of Bland- Hawthorn 
et al. (1991): these clearly have a higher resolution (and 
thus provide more apertures) but a more clumpy distribu- 
tion than the CO maps. We confirm the trend found by 
RW04 and clearly detect two different patterns, for points 
inside or outside 20". The scatter of the (V) versus (X) 
diagram (see Fig. 1121 is small for the outer region with a 
well-defined slope of 80 ± 2 km s" 1 kpc -1 for a PA of 90°, 
only slightly higher than the value found by RW04. How- 
ever, this slope strongly varies from 56 to 104 km s _1 kpc -1 
when we assume a PA varying from 84 to 95°, a depen- 
dency already emphasized by RW04 and Debattista (2003; 
see also Debattista & Williams 2004). The pattern speed 
of the inner region is rather badly constrained by our data, 
an expected result considering the presence of star forming 
regions and the influence of the ionization cone. Using the 
odd parts of (V) and {X) data points to minimize the scat- 
ter (Tremaine & Weinberg 1984), we find a lower limit value 
of 133 ± 12 km s _1 kpc -1 for PA= 90° (consistent with the 
value of RW04), and from 93 up to 178 km s" 1 kpc -1 for 
PAs again varying from 84 to 95°. 



5 N-BODY MODELLING 

We now turn to the dynamical modelling of the central re- 
gion of NGC 1068 using N-body and SPH simulations. This 
will be compared in Sect. 15.31 with the SAURON maps of the 
gas and stellar kinematics, as well as with the available HI, 
Ha and CO data (see Sect[3J. 

5.1 Methods 

We used PMSPH, the N-body code initially developed in 
Geneva (Pfenniger & Friedli 1993). The gravitational forces 
are computed with a particle-mesh method using a 3D polar 
grid with (Nr,N^,N z ) = (40,32,64) active cells. The hy- 
drodynamics equations are solved using the SPH technique 
(Friedli & Benz 1993). Since the radial spacing of the cells is 
logarithmic, the cell size increases from 15 pc at the centre 
to 383 pc at a radius of 2 kpc. 20 cells (i.e. half the number 
of radial cells) are used to described the region enclosed by 
SAURON field. The total grid has a radius of 100 kpc. All runs 
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Figure 12. Tremaine Weinberg method applied on the Ha Fabry 
Perot maps of Bland-Hawthorn et al. (1991). The averaged, den- 
sity weighted, velocity ( V) is plotted against the averaged, density 
weighted coordinates (X), along slits parallel to the line of nodes 
assuming PA= 90° . The best fit slopes for the outer points (open 
symbols, R > 20") and inner points (crosses, 8" < R < 15") are 
shown as a solid and dotted line, respectively. 

were performed using a total of 978 572 particles for the stel- 
lar part and 50000 for the gaseous one. The initial conditions 
(positions, velocities, velocity dispersions) were built from a 
Monte Carlo realization of a five-component axisymmetric 
model, as described below, and the total running times for 
each simulation varied from 500 to 2000 Myr. 

5.2 Mass model and initial conditions 

We constrained the initial axisymmetric mass model of 
NGC 1068 by using the available photometry and veloc- 
ity curve as constraints. We performed this in two steps. 
We first constrained the projected luminosity distribution 
of the galaxy by combining a high resolution near-infrared 
K band image (see Sect.|3J with the 2MASS K band image, 
the DSS I band image and a deep wide-field B band image 
(see Sect. HO)- The 

near-infrared images clearly reveal the 
bar within the central 20", and the optical images allowed 
us to estimate the contribution of the outer oval and disk 
structures (see Sect.[3J. We assumed the presence of a cen- 
tral spherical bulge and fitted a projected Plummer sphere 
to the high resolution NIR K band image. This component 
was subtracted from all images, which were then simply de- 
projected by assuming an inclination angle of i = 40°, and a 
two-dimensional distribution: this corresponds to a spatial 
scale factor of ~ 1.3 perpendicularly to the line of nodes. 
The deprojected images were found to be reasonably ap- 
proximated by the sum of four individual Miyamoto-Nagai 
components, each of them probing a different scale (see Ta- 
ble [TJ. 

We normalize each component, i.e. the four Miyamoto- 
Nagai, the Plummer sphere and one additional Miyamoto- 
Nagai component for the gas, so that the resulting circular 
velocity curve fits the observed HI velocity profile. This re- 
quires the overall mass-to-light ratio to increase by a factor 
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Figure 13. Circular velocity profiles of the Nbody+SPH simu- 
lations at t = (dotted lines) and at the end of the runs (solid 
lines) compared with observed velocity profiles, t = 850 for model 
A (thick lines) and t = 1440 for model B (thin lines). Symbols 
show the deprojected HI (filled squares), Ha (diamonds; from this 
paper) and CO (stars) velocity profiles (assuming an inclination 
of i = 40°). 
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of 10 between the central part and the outer most point of 
the HI curve (R = 180"). Since we know that the simulation 
will rapidly depart from its initial conditions, this procedure 
is only meant to produce a three dimensional mass distribu- 
tion which reproduces the radial mass gradient, the central 
bulge concentration and the large scale velocity curve rea- 
sonably well. We consider the contribution of a presumed 
central black hole to be part of the Plummer sphere. The 
initial radial profiles of the circular velocity are presented in 

Fig. ma 

We performed a number of runs (36), where we mainly 
varied the concentration of the different components, keep- 
ing the total mass (including the gas) roughly constant (see 
below). We chose to ignore the observed Ha velocity profile 
as a constrain for the initial conditions of our models as it 
exhibit strong non-circular motions (e.g., see Sect. 14.51 . We 
present only two models here (A & B), as they are typical 
of the runs performed and their initial circular velocity pro- 
files roughly bracket the observed HI profile of NGC 1068 
(Fig. 1131 . Both models have initial conditions which are 
consistent with a deprojected axisymmetric version of the 
available K band images. This is illustrated for model B in 
Fig. 1141 in which we show the projected mass profiles com- 
pared to 2MASS K band surface brightness profiles. The 
corresponding gravitational potential and density distribu- 
tions in the meridional plan are presented in the same figure. 



Figure 14. Initial conditions for model B. Top panels: poten- 
tial (left) and density (right) isocontours in the meridional plane 
(R, z). Contour steps are 0.1 and 0.5 in log respectively. Bottom 
panels: (projected) surface brightness profiles along the East- West 
(left) and South-North (right) direction extracted from the K 
band 2MASS images (open circles) and model B after projection 
using an inclination of 40° (solid lines). 



Table 1. Parameters of the initial conditions used in the N-body 
+ SPH simulations. M-N and P stand for Miyamoto-Nagai and 
Plummer components respectively. 



Component Type # Mass a b 

[10 9 M ] [kpc] [kpc] 



Model A 



Stars 


M-N 


1 


1.1 


0.255 


0.035 


Stars 


M-N 


2 


18.7 


1.0 


0.1 


Stars 


M-N 


3 


46.2 


5.0 


0.2 


Stars 


P 


4 


3.3 




0.1 


Gas 


M-N 




3.4 


6.0 


0.2 



Model B 



Stars 


M-N 


1 


2.1 


0.155 


0.035 


Stars 


M-N 


2 


18.7 


0.8 


0.2 


Stars 


M-N 


3 


46.2 


1.6 


0.4 


Stars 


P 


4 


3.3 




0.1 


Gas 


M-N 




3.4 


6.0 


0.2 



5.3 Models A and B 

A bar starts forming around 300 Myr, but becomes promi- 
nent only after 900 Myr for Model A. The bar forms much 
later on in Model B, this being mainly due to its significantly 
higher mass concentration. The overall evolution of Model 
B is much slower for the same reason. In simulations with 
even higher mass concentrations the disk was found to be 
stable against the formation of a bar. It is important to note 



that the central Plummer sphere has a characteristic scale 
(100 pc), 6 times the central resolution grid point (15 pc). 
After the bar formation, the gas wraps around in a low con- 
strast multi-arm spiral structure for a few 100 Myr. The gas 
then reacts quickly to the formation of the m — 2 perturba- 
tion and is radially redistributed. Part of the gaseous com- 
ponent flows towards the central 500 pc, exhibiting a time 
varying structure. We derived a snapshot of both simula- 
tions after the bar is well formed, the time of each snapshot 
being chosen (t = 850 for model A, and 1440 Myr for Model 
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B) such as to resemble the overall structure of the gaseous 
component in NGC 1068. We optimized the viewing angles 
for the bar to lie at a PA of 44.5° as in NGC 1068, and the 
zero velocity curve to be as close as possible to the observed 
one. This requires (for a fixed inclination of 40°) that the 
line of nodes is at a PA of ~ 83° and 84° for models A and 
B respectively. We have tested the dependence of this value 
on the time of the snapshot as well as on varying initial 
conditions: the resulting PA of the line of nodes varies then 
between 80 and 90°. 



5.4 Comparison with the observed stellar 
kinematics 

The simulations performed for this study are intended as 
a proof of concept, and not as perfect fits to the data. 
The kinematics are luminosity- weighted and we did not in- 
clude star formation or dust extinction in our simulations, 
although it is clear that young stars and dust significantly af- 
fect the surface brightness at the wavelengths of the SAUR0N 
observations. Furthermore, the bars formed in these simula- 
tions do not perfectly match the well-known near-infrared 
bar, and we could not reproduce the large-scale (North- 
South) oval in the photometry (the suggested signature of a 
large-scale bar). 

Model A fits the HI inner velocity profile better than 
model B (both models having similar circular velocities out- 
side 100"), but significantly fails to reproduce the first two 
stellar velocity moments (velocity and velocity dispersion), 
with values significantly below the observed data. We will 
therefore focus on the results from model B in the rest 
of the paper. The bar in model B has a pattern speed of 
~ 100 km s _1 kpc -1 , which leads to radii for the Inner Inner 
and Outer Inner Lindblad resonances (IILR, OILR), 4:1 ul- 
tra harmonic (UHR), Corotation (CR), and Outer Lindblad 
resonances (OLR) of 5'.'2, 10'.'5, 22'.'6, 32'.'2, 50'.'2, respec- 
tively (0.36, 0.73, 1.58, 2.25, 3.50 kpc respectively). This 
implies that the NIR bar ends well inside its CR, closer to 
its UHR. 

A comparison of the stellar kinematics of model B 
within the SAUR0N field of view is presented in Fig. 1151 it 
looks qualitatively similar to the observed SAUR0N maps. 
The twist in the stellar velocity field is less pronounced in 
the model, which nevertheless shows the clear signature of 
non-axisymmetry. A more quantative comparison is shown 
in Fig. 1161 where a few cuts at PA= 90° of the stellar veloc- 
ity fields are presented against the SAUR0N data. The overall 
agreement is quite good, with the exception of the outer 
parts of the central profile in which there is a discrepancy 
of about 20 km s" 1 : this is where we detected perturba- 
tions in the stellar kinematics fSect. l4~2"l . As mentioned al- 
ready, they follow the velocity differences between the H/3 
and [O ill] emission lines (Sect. l4~31l . The discrepancy be- 
tween the model and observed stellar kinematics therefore 
mostly reflects the presence of streaming motions within the 
outer spiral arms, where star formation is on-going. Model B 
additionally reproduces the morphology of the observed hs 
field, with the zero h$ curve almost following the major-axis 
of the bar, as in the SAUR0N map. There is only a slight cen- 
tral depression in the hi map of model B, less pronounced 
than in the data. The model does also not reproduce the 



200 



-200 




x [arcaec] 



Figure 16. Comparison between stellar kinematic profiles ex- 
tracted from the two-dimensional SAUR0N velocity field (solid 
lines) and from model B at i = 1440 Myr (dashed lines). Four 
cuts along PA= 90° are shown, with offsets, from top to bottom, 
of 10, 5, 0, -5, and -10 arcseconds respectively. 



ring-like enhancement of hi observed both in the SAUR0N 
and long-slit data. 

We can now compare the obtained stellar kinematics 
with the more extended long-slit kinematics obtained from 
the data of Sh+03: a comparison is presented in Fig. 1171 The 
agreement is again very reasonable, including the higher or- 
der Gauss-Hermite moments. The kinematic profiles of the 
simulations do not however exhibit the characteristic dis- 
persion peaks at a radius of ~ 22" on either sides of the 
center along the major-axis. These peaks correspond to a 
transition region at the end of the bar, where the slope of 
the velocity gradient also changes, and are close to the UHR 
in model B. The simulations also lack the central plateau or 
depression in the dispersion profile observed in the central 
5" of NGC 1068. 



6 DISCUSSION 

6.1 The location of resonances 

Both the molecular and ionized gas kinematics exhibit 
strong departures from circular motion (see e.g., Fig. 1111 . 
and the presence of the near- infrared bar in the inner 1.5 kpc 
is undisputable. The Ha (Fig. 1 1 UP and CO velocity fields 
(S+00) provide us with an important clue regarding the 
kinematics of the spiral arms: the isovelocities clearly bend 
outwards which implies an inward streaming motion (e.g. 
see the spiral arm velocity contours at about 10" South, 5" 
West). A similar argument was initially emphasized by Yuan 
& Kuo (1998), although they reached the opposite conclu- 
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Figure 15. Luminosity distribution (top left panel, assuming a constant mass-to-light ratio), and observed stellar kinematics from the 
model B at t = 1440 (using an inclination of i = 40°). The panels correspond to the ones in Fig.|S] Boxes in the top right corner of the 
panels show the ranges covered by the color bars. 



sion as they relied at the time on the CO data of Heifer & 
Blitz (1995) which lacked spatial resolution and seemed to 
show an inward bending 3 . The inward streaming therefore 
requires the gas to be located inside the corotation (CR) or 
outside the outer Lindblad resonance (OLR) of the corre- 
sponding density wave, a conclusion also reached by S+00 
and RW04 (and opposite to the one of Yuan & Kuo 1998, 
but see above). 

The presence of a large-scale tumbling component with 
a radius of about 9 kpc is debatable, although a very signif- 
icant elongation is clearly observed in the photometry at a 
PA near the minor-axis, inconsistent with an axisymmetric 
outer disk. S+00 first estimated the pattern speed of the 
outer oval 4 to be ft p ~ 35 km s -1 kpc -1 . Assuming the 
ILR of the primary bar corresponds to the CR of the in- 
ner (secondary) bar, they determine its pattern speed to be 



3 Yuan & Kuo (1998) also noted that they "would not be sur- 
prised" to get a different result when a more reliable rotation 
curve becomes available. 

4 S+00 quotes a value of O p = 20 km s — 1 kpc -1 in their paper, 
although, as mentioned by RW04, they actually used a value of 
35 km s — 1 kpc -1 . 



£l s ~ 140 km s -1 kpc -1 . This would imply that the gaseous 
spiral arms lie between the CR and ILR of the primary bar. 
These arms would however lie outside the CR of the sec- 
ondary bar, and the assumed value for fi s would make the 
near-infrared bar extend as far as its CR. Recent simulations 
of double bars (Rautiainen, Salo & Laurikaine 2004) seem 
indeed to favor CR/ILR coupling (but see Moellenhoff et 
al. 1995), but with short secondary bars, ending inside the 
4:1 Ultra Harmonic Resonance (UHR). This requires the 
pattern speeds for both the primary and secondary bar to 
be lower. 

Assuming the PA of the line of nodes is between 84 
and 90° (Sect. 15. 41 . a lower limit for the pattern speed 
for the inner bar is estimated to be between 93 and 
133 km s -1 kpc -1 (Tremaine- Weinberg method on the Ha 
data, see Sect. 14.611 . while the bar that forms in Model B 
has a pattern speed of 100 km s -1 kpc - (Sect. l5~4t . These 
values would be consistent with the inner bar ending inside 
its UHR as discussed above. If we believe the outer oval cor- 
responds to a third tumbling structure (S+00), a CR/ILR 
coupling between the two bars would then imply a primary 
bar with Q p ~ 25 km s _1 kpc -1 , lower than but still consis- 
tent with the value advocated by S+00. 
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Figure 17. Comparison between stellar kinematic profiles ex- 
tracted from the Sh+03 (crosses) and from model B at t = 
1440 Myr (solid lines). From top to bottom: mean velocity, ve- 
locity dispersion, hg and /14. The model kinematics have been 
averaged over a slit width of 3". 

As mentioned above, the N body simulations conducted 
in this paper lead to a PA for the lines of nodes in the range 
83-90° with a preferred value of 84° . This is consistent with 
the photometric PA of the outer disk, and with the out- 
put of the tilted ring gas modelling (Sect. 14. 51 . It is not 
consistent with the HI kinematic axis determined by Brinks 
(1997). Brinks however noticed a gradual change in the po- 
sition angle of the HI kinematic major axis, which he inter- 
preted as the sign of a mild warping. The HI map exhibits 
a weak spiral arm or ring-like structure in the radius range 
130-180" the density decreasing then inwards (down to the 
inner HI ring present between 30" and 80"; Brinks 1997). 
If the outer oval is, as mentioned above, a weak bar tum- 
bling at 25 km s _1 kpc -1 , its Outer Lindblad Resonance 
would be located at a radius of 140" (10 kpc), within the 
outer HI structure. At the OLR, we expect periodic orbits to 
be elongated perpendicularly to the major-axis of the bar. 
The observed radial variation in the HI kinematics could 
therefore be partly explained by the resulting non-circular 
motions. This should however be confirmed by a detailed 
kinematic HI study in the outer region of NGC 1068. 

We finally turn to the derivation of the pattern speed 
associated with the outer spiral arms. Our estimate us- 
ing the Tremaine- Weinberg technique is between 56 and 
80 km s -1 kpc -1 (for PAs between 84 and 90°), consistent 
with the values given by RW04. This is however significantly 
lower than the estimate obtained via the harmonic decom- 
position (Sect. 1131 . which provides a value close to the esti- 



mate of the pattern speed of the NIR bar. If the outer spiral 
and the inner bar share the same pattern speed, it would 
hint at the bar being the main driver. But the outer spiral 
and the inner bar having different pattern speeds would not 
be unexpected as modes can couple non-linearly (Masset & 
Tagger 1997). Such a process is witnessed in some simula- 
tions of double bars (see e.g., Rautiainen et al. 2004) which 
exhibit up to four different m — 2 modes. It would however 
not be wise, considering the large uncertainties in the deter- 
mination of the tumbling frequencies, to speculate further if 
non-linear coupling of modes are present in NGC 1068. The 
pattern speed of the outer spiral is in all cases much higher 
than the one of the presumed large-scale bar. 

6.2 Signatures of the bar 

Bureau & Athanassoula (2005) drew attention to a num- 
ber of the stellar kinematics signatures of a bar using pure 
N body simulations: a double hump velocity profile, corre- 
lated stellar V and fi3 profiles within the bar, secondary 
peaks for the dispersion a near the end of the bar, a flat 
central minimum in the /14 profile followed by a significant 
rise and a decline at larger radii. These features were how- 
ever obtained for rather high inclination (i > 75°), and for 
a galaxy with a very prominent outer stellar ring. There 
are secondary peaks in the observed dispersion profiles of 
NGC 1068 at a radius of 22" for a PA of 80° (see e.g., Fig. lTTjl . 
and a flat minimum in hi which corresponds to the flat cen- 
tral part in a (within a radius of 5"). But, although there 
is a local minimum in V corresponding to the secondary 
a maxima, the double hump nature of the velocity profile 
is not very pronounced: the outer stellar velocity curve is 
nearly flat with a level comparable to the inner maximum 
of ~ 120 km s _1 at radii between 10 and 15". We further- 
more observe /13 profiles which are clearly anti-correlated 
with V within the extent of the bar: this means that we de- 
tect a tail of low velocity stars superimposed on the rapidly 
rotating ones along the line of nodes. 

The high resolution K band image of NGC 1068 ob- 
tained by Peletier et al. (1999) shows the inner bar with 
rather pinched ends at a radius of ~ 16". These can be inter- 
preted as the tips of the xi orbits sustaining the bar. There is 
a clear plateau in the surface brightness profile around that 
radius, which is also the witness of a radial (and possible 
vertical) redistribution of the stellar component (Combes et 
al. 1990). The fact that we do not reproduce the higher stel- 
lar dispersion (and lower stellar velocity) around a radius 
of 22" (PA= 80° , Fig. E3 with our simulations may be due 
to the lack of star formation in our models. This region in 
NGC 1068 is the site of intense star formation, and corre- 
sponds to the UHR of model B: stars often form near the 
UHR (Schwarz 1984, Buta & Combes 1996), and gathering 
new stars with orbits at the UHR (4:1 resonance) may be 
the cause of an increase in the stellar velocity dispersion. 

6.3 Driving gas towards the central 300 pc 

Star formation could also be the cause of the discrepancy be- 
tween the rather flat (or depressed) central dispersion profile 
of NGC 1068 and that of model B. The inner bar is able to 
drive gas inwards, as illustrated by some numerical simu- 
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lations of double barred galaxies (see Model III of Rauti- 
ainen et al. 2004) although this is certainly not a generic 
result for double bars (e.g., Maciejewski et al. 2002). If we 
assume that the NIR bar ends well inside its CR, it then im- 
plies the presence of an Inner Lindblad Resonance around 
5". This would roughly coincide with the extent of the stel- 
lar velocity decoupling and cr-drop, and support the generic 
scenario outlined by Emsellem et al. (2001) and Wozniak et 
al. (2003): the inner NIR bar has driven gas within its Inner 
Lindblad Resonance, thus forming a relatively cold stellar 
system within the central few hundred parcsec. 

We are indeed witnessing radial gas inflow in the central 
kpc in our numerical simulations. The inflow rate in model B 
ranges from 10 -2 to a few MQ/yr within the region of the 
bar, with an average of about 5 10 -2 Mo/yr. This gaseous 
mass inflow increases significantly as the bar gets stronger. 
The net mass increase driven by the bar 300 Myr after its 
formation is ~ 5 10 7 Mq inside 5" (350 pc), which represents 
about 1% of the total enclosed mass. All these values should 
be taken as upper limits as e.g., the simulations do not in- 
clude star formation. In order to make the link with still 
smaller scales, high resolution two-dimensional mapping of 
the stellar kinematics and populations in the inner few arc- 
seconds would be required - a considerable challenge given 
the active and disturbed nuclear environment of NGC 1068. 



NGC 1068. We focus on one model which has a gravita- 
tional potential consistent with the HI velocity profile, and 
in which a bar forms. After projection with the relevant 
viewing angles, the resulting stellar kinematics successfully 
reproduce a number of properties observed in the SAUR0N 
and long-slit kinematics, including the /13 and /14 profiles. 
There are however some discrepancies between the model 
and the observations, which could be partly explained by 
the exclusion of star formation and the lack of a primary 
large-scale bar in our simulations. 

We finally briefly discuss the possibility of gas fueling 
within the inner bar. The numerical simulations suggest that 
the inner bar could drive a significant amount of gas down 
to a scale of ~ 300 pc (its ILR). This would be consistent 
with the interpretation of the cr-drop (Emsellem et al. 2001) 
in NGC 1068 being the result of central gas accretion fol- 
lowed by an episode of star formation (Emsellem et al. 2001; 
Wozniak et al. 2003) . We should however wait for a detailed 
view at the stellar kinematics and populations within the 
inner 300 pc to conclude on this issue, e.g. with an integral- 
field spectrograph in the near-infrared such as SINFONI at 
the VLT. It would also be important to examine a sample 
of carefully chosen galaxies, as to allow a statistical anal- 
ysis and to assess the relative contributions of the various 
physical mechanisms at play. 



7 CONCLUSIONS 

We have successfully obtained stellar and gaseous kinemat- 
ics for NGC 1068 over a large two-dimensional field covering 
the entire near-infrared inner (secondary) bar. Although our 
spectra exhibit numerous emission and absorption features, 
we have been able to disentangle the two components, and 
derive the kinematics independently. Our SAUR0N data re- 
veal a regular stellar velocity field that rules out an offset 
stellar system as claimed by Garcia-Lorenzo et al. (1998). 
The SAURON two-dimensional stellar kinematics exhibit the 
signatures of the presence of the inner bar: a twist in the ve- 
locity field, and an /13 field elongated along the bar. These 
features are also retrieved in a qualitative comparison with 
N-body + SPH simulations. 

We detect a kinematic decoupling of the central 350 pc 
(5") in a reanalysis of previously published long-slit data 
(Sh+03). We first observe a flattening of the h?, profile indi- 
cating a change of orbital structure within this region. We 
then reveal a flattening or a potential drop in the central 
part of the stellar velocity dispersion profile. The kinemat- 
ics of the ionized H/3 and [O ill] gas both show strong in- 
wards streaming motions. Differences up to 100 km s _1 are 
observed in the ionised gas kinematics inferred from the H/3 
and [O ill] lines. The CO and Ha distributions are not coin- 
cident, a difference which can be explained by the effect of 
the on-going star formation, and dust extinction. 

We confirm the presence of two different pattern speeds 
for the region inside and outside of the inner bar (RW04), 
applying the Tremaine- Weinberg method to the Ha veloc- 
ity field. The pattern speed of the inner bar is however only 
weakly constrained. Note that the outer oval could corre- 
spond to a third tumbling structure, as mentioned by S+00. 

We then performed numerical simulations with initial 
axisymmetric conditions approximating the photometry of 
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